Optimized suppression of coherent noise from seismic data using the Karhunen-Loeve 

transform 
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Signals obtained in land seismic surveys are usually contaminated with coherent noise, among 
which the ground roll (Rayleigh surface waves) is of major concern for it can severely degrade the 
quality of the information obtained from the seismic record. Properly suppressing the ground roll 
from seismic data is not only of great practical importance but also remains a scientific challenge. 
Here we propose an optimized filter based on the Karhunen-Loeve transform for processing seismic 
data contaminated with ground roll. In our method, the contaminated region of the seismic record, 
to be processed by the filter, is selected in such way so as to correspond to the maximum of a 
properly defined coherence index. The main advantages of the method are that the ground roll is 
suppressed with negligible distortion of the remanent reflection signals and that the filtering can 
be performed on the computer in a largely unsupervised manner. The method has been devised to 
filter seismic data, however it could also be relevant for other applications where localized coherent 
structures, embedded in a complex spatiotemporal dynamics, need to be identified in a more refined 
way. 

PACS numbers: 93.85.+q,91.30.Dk ,43.60.Wy,43.60.Cg 
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I. INTRODUCTION 

Locating oil reservoirs that are economically viable is 
one of the main problems in the petroleum industry. This 
task is primarily undertaken through seismic exploration, 
where explosive sources generate seismic waves whose re- 
flections at the different geological layers are recorded at 
the ground or sea level by acoustic sensors (geophones 
or hydrophones). These seismic signals, which are later 
processed to reveal information about possible oil oc- 
currences, are often contaminated by noise and properly 
cleaning the data is therefore of paramount importance 
. In particular, the design of efficient filters to suppress 
noise that shows coherence in space and time (and often 
appears stronger in magnitude than the desired signal) 
remains a scientific challenge for which novel concepts 
and methods are required. In addition, the filtering tools 
developed to treat such kind of noise may also find rel- 
evant applications in other physical problems where co- 
herent structures evolving in a complex spatiotemporal 
dynamics need to identified properly. 

In land seismic surveys, the seismic sources generate 
various type of surface waves which are regarded as noise 
since they do not contain information from the deeper 
subsurface. This so-called coherent noise represents a se- 
rious hurdle in the processing of the seismic data since 
it may overwhelm the reflection signal, thus severely de- 
grading the quality of the information that can be ob- 
tained from the data. A source-generated noise of partic- 
ular concern is the ground roll, which is the main type of 
coherent noise in land seismic records and is commonly 
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much stronger in amplitude than the reflected signals. 
Ground roll are surface waves whose vertical components 
are Rayleigh-type dispersive waves, with low frequency 
and low phase and group velocities. 

An example of seismic data contaminated by ground 
roll is shown in Fig. ^ This seismic section consists of 
land-based data with 96 traces (one for each geophone) 
and 1001 samples per trace. A typical trace is shown in 
Fig. |21 corresponding to geophone 58. The image shown 
in Fig. ^ was created from the 96 traces using a stan- 
dard imaging technique. The horizontal axis in this figure 
corresponds to the offset distance between source and re- 
ceiver and the vertical axis represents time, with the ori- 
gin located at the upper-left corner. The maximum offset 
is 475 m (the distance between geophones being 5 m) and 
the maximum time is 1000 ms. The gray levels in Fig. 1 
change linearly from black to white as the amplitude of 
the seismic signal varies from minimum to maximum. 
Owing to its dispersive nature, the ground roll appears 
in a seismic image as a characteristic fan-like structure, 
which is clearly visible in Fig. ^ The data shown in this 
figure was provided by the Brazilian Petroleum Company 
(PETROBRAS). 

Standard methods for suppressing ground roll include 
one-dimensional high-pass filtering and two-dimensional 
f-k filtering [|. Such "global" filters are based on the 
elimination of specific frequencies and have the disad- 
vantage that they also affect the uncontaminated part 
of the signal. Recently, "local" filters for suppressing 
the ground roll have been proposed using the Karhunen- 
Loeve transform 0, 0] and the wavelet transform 0, . 
The Wiener-Levinson algorithm has also been applied to 
extract the ground roll [6j . 

Filters based on the Karhunen-Loeve (KL) transform 
are particularly interesting because of the adaptativity 
of the KL expansion, meaning that the original signal 
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FIG. 1: A space-time plot of seismic data. The horizontal axis 
represents the offset distance and the vertical axis indicates 
time. The origin is at the upper-left corner, and the maximum 
offset and time are 475 m and 1000 ms, respectively. The gray 
scale is such that black (white) corresponds to the minimum 
(maximum) amplitude of the seismic signal. The ground roll 
noise appears as downward oblique lines. 




250 500 750 1000 

time 



FIG. 2: Seismic signal recorded by a single geophone (trace 
58). The amplitude is in arbitrary units and time in ms. 

is decomposed in a basis that is obtained directly from 
the empirical data, unlike Fourier and wavelet transforms 
which use prescribed basis functions. The KL transform 
is a mathematical procedure (also known as proper or- 
thogonal decomposition, empirical orthogonal function 
decomposition, principal component analysis, and singu- 
lar value decomposition) whereby any complicated data 
set can be optimally decomposed into a finite, and often 
small, number of modes (called proper orthogonal modes, 
empirical orthogonal functions, principal components or 
eigenimages) which are obtained from the eigenvectors 



of the data autocorrelation matrix. In applying the KL 
transform to suppress the ground roll, one must first map 
the contaminated region of the seismic record into a hor- 
izontal rectangular region. This transformed region is 
then decomposed with the KL transform and the first 
few principal components are removed to extract the co- 
herent noise, after which the filtered data is inversely 
mapped back into the original seismic section. The ad- 
vantage of this method is that the noise is suppressed 
with negligible distortion of the reflection signals, for only 
the data within the selected region is actually processed 
by the filter. Earlier versions of the KL filter 0, have 
however one serious drawback, namely, the fact that the 
region to be filtered must be picked by hand — a proce- 
dure that not only can be labor intensive but also relies 
on good judgment of the person performing the filtering. 

In this article we propose a significant improvement of 
the KL filtering method, in which the region to be filtered 
is selected automatically as an optimization procedure. 
We introduce a novel quantity, namely, the coherence in- 
dex CI, which gives a measure of the amount of energy 
contained in the most coherent modes for a given selected 
region. The optimal region is then chosen as that that 
gives the maximum CI. We emphasize that introduc- 
ing a quantitative criterion for selecting the 'best' region 
to be filtered has the considerable advantage of yielding 
a largely unsupervised scheme for demarcating and effi- 
ciently suppressing the ground roll. 

Although our main motivation here concerns the sup- 
pression of coherent noise in seismic data, we should like 
to remark that our method may be applicable to other 
problems where coherent structures embedded in a com- 
plex spatiotemporal dynamics need to be identified or 
characterized in a more refined way. For example, the KL 
transform has been recently used to identify and extract 
spatial features from a complex spatiotemporal evolution 
in combustion experiment 0, H, @ ■ A related method — 
the so-called biorthogonal decomposition — has also been 
applied to characterize spatiotemporal chaos and iden- 
tified structures [l(J, Ej as weu as identify changes in 
the dynamical complexity, and the spatial coherence of 
a multimode laser [l^l- We thus envision that our opti- 
mized KL filter may find applications in these and related 
problems of coherent structures in complex spatiotempo- 
ral dynamics. 

The article is organized as follows. In Sec. [D] wc define 
the Karhunen-Loeve transform, describe its main proper- 
ties, and discuss its relation to the singular value decom- 
position of matrices. In Sec. IHII we present the KL filter 
and a novel optimization procedure to select the noise- 
contaminated region to be parsed through the filter. The 
results of our optimized filter when applied to the data 
shown in Fig. ^ are presented in Sec. IIVI Our main con- 
clusions are summarized in Sec. In Appendixes A 
and B we briefly discuss, for completeness, the relation 
between the KL transform and two other similar pro- 
cedures known as proper orthogonal decomposition (or 
empirical orthogonal function expansion) and principal 
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component analysis. 



II. THE KARHUNEN-LOEVE TRANSFORM 

A. Definition and main properties 

Consider a multichannel seismic data consisting of m 
traces with n samples per trace represented by a m x n 
matrix A, so that the element Aij of the data matrix cor- 
responds to the amplitude registered at the ith geophonc 
at time j. For dcfinitcness, let us assume that m < n, as 
is usually the case. We also assume for simplicity that 
the matrix A has full rank, i.e., r = m, where r denotes 
the rank of A. Letting the vectors Xi and y~j denote the 
elements of the ith row and the jth column of A, respec- 
tively, we can write 



A = (yi y 2 •■■ y n ) 



•?2 



(1) 



us denote by the Xii * — 1, m, the elements of the ith 
row of the KL matrix 'J, that is, 



X2 



(7) 



V Xr, 



Then JBJ can be written as 



i=l 



(8) 



where it is implied matrix multiplication between the col- 
umn vector vii and the row vector The eigenvectors 
Ui are called empirical eigenvectors, proper orthogonal 
modes, or KL modes. 

As discussed in Appendix [S] the total energy E of the 
data can be defined as the sum of all eigenvalues, 



E : 



(9) 



With the above notation we have 



VH 



(2) 



where ay denotes the jth element of the vector Oj . (To 
avoid risk of confusion matrix elements will always be 
denoted by capital letters, so that a small-cap symbol 
with two subscripts indicates vector elements.) 

Next consider the following m x m symmetric matrix 



AA 



(3) 



where the superscript t denotes matrix transposition. It 
is a well known fact from linear algebra that matrices of 
the form pfll. also called covariance matrices, are positive 
definite . Let us then arrange the eigenvalues A; of T 
in non-ascending order, i.e., Ai > A2 > ... > A m > 0, and 
let Ui be the corresponding (normalized) eigenvectors. 

The Karhunen-Loeve (KL) transform of the data ma- 
trix A is defined as the m x n matrix ^ given by 



f = WA, 



(4) 



where the columns of the matrix U are the eigenvectors 
of T: 



U = (ui u 2 



(5) 



The original data can be recovered from the KL trans- 
form ^ by the inverse relation 



A = [/*. 



(6) 



We refer to this equation as the KL expansion of the data 
matrix A. To render such an expansion more explicit let 



so that Ai can be interpreted as the energy captured by 
the ith empirical eigenvector u,. We thus define the rel- 
ative energy Ei in the ith KL mode as 



E; 



(10) 



We note furthermore that since T is a covariance-like ma- 
trix its eigenvalues A.; can also be interpreted as the vari- 
ance of the respective principal component Uf, see Ap- 
pendix[H]for more details on this interpretation. We thus 
say that the higher Ai the more coherent the KL mode 
Ui is. In this context, the most energetic modes are iden- 
tified with the most coherent ones and vice-versa. 

An important property of the KL expansion is that it is 
'optimal' in the following sense: if we form the matrix 
by keeping the first k rows of \& and setting the remaining 
to — k rows to zero, then the matrix Ak given by 



A k = UVi 



(11) 



is the best approximation to A by a matrix of rank k < m 
in the Frobenius norm (the square root of the sum of 
the squares of all matrix elements) fl3l |. This optimal- 
ly property of the KL expansion lies at the heart of its 
applications in data compression 01 an d dimensionality 
reduction [13j, for it allows to approximate the original 
data A by a smaller matrix Ak with minimum loss of in- 
formation (in the above sense). Another interpretation 
of relation l|llf) is that it gives a low-lass filter [15( , for in 
this case only the first k KL modes are retained in the 
filtered data Ak- 

On the other hand, if the relevant signal in the appli- 
cation at hand is contaminated with coherent noise, as 
is the case of the ground roll in seismic data, one can 
use the KL transform to remove efficiently such noise by 
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constructing a high-pass filter. Indeed, if we form the 
matrix ^f' k by setting to zero the first k rows of W and 
keeping the remaining ones, then the matrix A'j~ given 
by 

A' k = U% (12) 

is a filtered version of A where the first k 'most coher- 
ent' modes have been removed. However, if the noise is 
localized in space and time it is best to apply the filter 
only to the contaminated part of the signal. In previous 
versions of the KL filter the choice of the region to be 
parsed through the filter was made a priori, according to 
the best judgement of the person carrying out the filter- 
ing, thus lending a considerable degree of subjectivity to 
the process. In the next section, we will show how one 
can use the KL expansion to implement an automated 
filter where the undesirable coherent structure can be 
'optimally' identified and removed. 

Before going into that, however, we shall briefly discuss 
below an important connection between the KL trans- 
form and an analogous mathematical procedure known 
as the singular value decomposition of matrices. Read- 
ers already knowledgeable about the equivalence between 
these two formalisms (or more interested in the specific 
application of the KL transform to filter coherent noise) 
may skip the remainder of this section without loss of 
continuity. 

B. Relation to Singular Value Decomposition 

We recall that the singular value decomposition (SVD) 
of any m x n matrix A, with m < n, is given by the 
following expression: 

A = UY,V\ (13) 

where U is as defined in J5J|, Eisamxm diagonal matrix 
with elements <7j = vXj , the so-called singular values of 
A, and V is a m x n matrix whose columns correspond to 
the m eigenvectors {tfj} of the matrix A 1 A with nonzero 
eigenvalues. The SVD allows us to rewrite the matrix A 
as a sum of matrices of unitary rank: 

tci m 

A = Y1 (Ti( 5 i = E ( 14 ) 
»=i i=i 

In the context of i mag e processing the matrices Qi are 
called eigenimages jig. 

Now, comparing JHJ with (|13fl we see that the KL 
transform ^> is related to the SVD matrices £ and V 
by the following relation 

* = Sy*, (15) 

so that the row vectors Xi of ^ are given in terms of the 
singular values cr* and the vectors Vi by 

Xi = o- l v i t . (16) 
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FIG. 3: Schematic diagram for demarcating the ground roll 
on a seismic section and the corresponding rectangular sectors 
obtained by applying a linear map. 

It thus follows that the decomposition in eigenimages 
seen in (f 141) is precisely the KL expansion given in 
Furthermore the approximation A^ given in (|ll|l can be 
written in terms of eigenimages as 

k 

A k = J2°iQi- ( 17 ) 

i=l 

Similarly, the filtered data A' k shown in 112(1 reads in 
terms of eigenimages: 

m 

4= E ( 18 ) 

i=k+l 

The SVD provides an efficient way to compute the KL 
transform, and we shall use this method in the numerical 
procedures described in the paper. 

III. THE OPTIMIZED KL FILTER 

As already mentioned, owing to its dispersive nature 
the ground-roll noise appears in a seismic image as a 
typical fan-like coherent structure. This space-time lo- 
calization of the ground roll allows us to apply a sort of 
'surgical procedure' to suppress the noise, leaving intact 
the uncontaminated region. To do that, we first pick lines 
to demarcate the start and end of the ground roll and, if 
necessary, intermediate lines to demarcate different wave- 
trains, as indicated schematically in Fig. [3] In this figure 
we have for simplicity used straight lines to demarcate 
the sectors but more general alignment functions, such 
as segmented straight lines, can also be chosen 0,0- To 
make our discussion as general as possible, let us assume 
that we have a set of N parameters i = 1, N, de- 
scribing our alignment functions. For instance, in Fig. |3| 
the parameters {9i} would correspond to the coefficients 
of the straight lines defining each sector. 
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Once the region contaminated by the ground roll has 
been demarcated, we map each sector onto a horizon- 
tal rectangular region by shifting and stretching along 
the time axis; see Fig. [3] The data points between the 
top and bottom lines in each sector is mapped into the 
corresponding new rectangular domain, with the map- 
ping being carried out via a cubic convolution interpola- 
tion technique |17| . After this alignment procedure the 
ground roll events will become approximately horizontal, 
favoring its decomposition in a smaller space. Since any 
given transformed sector has a rectangular shape it can 
be represented by a matrix, which in turn can be decom- 
posed in empirical orthogonal modes (eigenimages) using 
the KL transform. The first few modes, which contain 
most of the ground roll, are then subtracted to extract 
the coherent noise. The resulting data for each trans- 
formed sector is finally subjected to the corresponding 
inverse mapping to compensate for the original forward 
mapping. This leaves the uncontaminated data (lying 
outside the demarcated sectors) unaffected by the whole 
filtering procedure. 

The KL filter described above has indeed shown good 
performance in suppressing source-generated noise from 
seismic data 0, 13 • The method has however the draw- 
back that the region to be filtered must be picked by 
hand, which renders the analysis somewhat subjective. 
In order to overcome this difficulty, it would be desirable 
to have a quantitative criterion based on which one could 
decide what is the 'best choice' for the parameters {Oi] 
describing the alignment functions. In what follows, we 
propose an optimization procedure whereby the region to 
be filtered can be selected automatically, once the generic 
form of the alignment functions is prescribed. 

Suppose we have chosen I sectors to demarcate the dif- 
ferent wavetrains in the contaminated region of the orig- 
inal data, and let {9\, ...,6n} be the set of parameters 
characterizing the respective alignment functions that de- 
fine these sectors. Let us denote by Ak, k = 1, the 
matrix representing the kth transformed sector obtained 
from the linear mapping of the respective original sec- 
tor, as discussed above. For each transformed sector Ak 
we then compute its KL transform and calculate the co- 
herence index Clk for this sector, defined as the relative 
energy contained in its first KL mode: 

Ch = = (19) 

where Af are the eigenvalues of the correlation matrix 
Tfc = AkA l k and Tk is the rank of A^. Such as defined 
above, Clk represents the relative weight of the most 
coherent mode in the KL expansion of the transformed 
sector Ak- (A quantity analogous to our CI is known in 
the oceanography literature as the similarity index |18|.') 

Next we introduce an overall coherence index 
CI(6±, ...,6m) for the entire demarcated region, defined 



as the average coherence index of all sectors: 

1 1 

CI(6 1 ,...,8 N ) = 1 J2ci k . (20) 

k=l 

As the name suggests, the coherence index CI is a mea- 
sure of the amount of 'coherent energy' contained in 
the chosen demarcated region given by the parameters 
{8i}iLi- Thus, the higher CI the larger the energy con- 
tained in the most coherent modes in that region. For the 
purpose of filtering coherent noise it is therefore mostly 
favorable to pick the region with the largest possible CI. 
We thus propose the following criterion to select the op- 
timal region to be filtered: vary the parameters {9{\ over 
some appropriate range and then choose the values 9* 
that maximize the coherence index CI, that is, 

CI(9l, ...,9* N ) ^ max[CI(9 1 , ...,9 N )} . (21) 

Once we have selected the optimal region, given by the 
parameters {9*}fL 1 , we then simply apply the KL fil- 
ter to this region as already discussed: we remove the 
first few eigenimages from each transformed sector and 
inversely map the data back into the original sectors, so 
as to obtain the final filtered image. In the next section 
we will apply our optimized KL filtering procedure to the 
seismic data shown in Fig. ^ 

IV. RESULTS 

Here we illustrate how our optimized KL filter works 
by applying it to the seismic data shown in Fig. ^ In 
this case, it suffices to choose only one sector to demar- 
cate the entire region contaminated by the ground roll. 
This means that we have to prescribe only two alignment 
functions, corresponding to the uppermost and lower- 
most straight lines (lines AB and CD, respectively) in 
Fig. |21 To reduce further the number of free parame- 
ters in the problem, let us keep the leftmost point of 
the upper line (point A in Fig. [3} fixed to the origin, so 
that the coordinates {ia^a) of point A are set to (0,0), 
while allowing the point B to move freely up or down 
within certain range; see below. Similarly, we shall keep 
the rightmost point of the lower line (point C in Fig. |3J) 
pinned at a point {ic,jc)i where ic — 95 and jc is cho- 
sen so that the entire ground roll wavetrain is above this 
point. The other endpoint of the lower demarcation line 
(point D in Fig. [3J| is allowed to vary freely. With such 
restrictions, we are left with only two free parameters, 
namely, the angles 9\ and 9^ that the upper and lower 
demarcation lines make with the horizontal axis. So re- 
ducing the dimensionality of our parameter space allows 
us to visualize the coherence index CI[9\,92) as a 2D 
surface. For the case in hand, it is more convenient how- 
ever to express CI not as a function of the angles 9\ and 
#2 but in terms of two other new parameters introduced 
below. 



Let the coordinates of point B, which defines the right 
endpoint of the upper demarcation line in Fig. [21 be given 
by (13,33)1 where is = 95. In our optimization proce- 
dure we let point B move along the right edge of the seis- 
mic section by allowing the coordinate js to vary from 
a minimum value jB min to a maximum value js max , so 
that we can write 

3b = j Bmin + kA B , k = 0,l, ...,N B (22) 

where Nb is the number of intermediate sampling 
points between ]B min and js mox , and Ab = (jB mam — 
3B min ) /Nb- Similarly, for the coordinates (in, jo) of 
point D in Fig. |3 which is the moving endpoint of the 
lower straight line, we have iu = and 

]D=j Dmm +lA D , l = 0,l,...,N D , (23) 

where is the number of sampling points between 
3D min and jD mam , and A D = (jc mQX - 3D min )/N D . 




FIG. 4: The coherence index C7 as a function of the indices 
k and I that define the demarcation lines; see text. 

For each choice of k and I in <(2*2^l and l(2*3"|l , we apply the 
procedure described in the previous section and obtain 
the coherence index CI(k, I) of the corresponding region. 
In Fig. 0] we show the energy surface CI(k,l), for the 
case in which j's TOin = 280, jB max = 600, jc = 864, 
3D min = 0, j Dnax = 576, and Ng — Nb — 64. We 
see in this figure that CI possesses a sharp peak, thus 
showing that this criterion is indeed quite discriminating 
with respect to the positioning of the lines demarcating 
the region contaminated by the ground roll. The global 
maximum of CI in Fig.0|is located at k = 42 and I = 24, 
and in Fig. [5Ji we show the transformed sector obtained 
from the linear mapping of this optimal region. In this 
figure one clearly sees that the ground roll wavetrains 
appear mostly as horizontal events. In Fig. [3Jd we present 
the first eigenimage of the data shown in Fig. which 
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FIG. 5: a) The selected region in the new domain; b) its first 
eigenimage; c) the second eigenimage; and d) the result after 
subtracting the first eigenimage. 



corresponds to about 33% of the total energy of the image 
in Fig. Ek, as can be seen in Fig. where we plot the 
relative energy Ei captured by the first 10 eigenimages. 
The second eigenimage, shown in Fig. El 3 , captures about 
10% of the total energy, with each successively higher 
mode contributing successively less to the total energy; 
see Fig. In Fig. [S|I we give the result of removing 
the first KL mode (Fig. [Sfc>) from Fig. [SJu It is clear 
in Fig. 0i that by removing only the first eigenimage the 
main horizontal events (corresponding to the ground roll) 
have already been greatly suppressed. 

Performing the inverse mapping of the image shown in 
Fig. yields the data seen in the region between the two 
white lines in Fig.[7Ji, which shows the final filtered image 
for this case (i.e., after removing the first KL mode from 
the transformed region). We see that the ground roll in- 
side the demarcated region in Fig. has been consider- 
ably suppressed, while the uncontaminated signal (lying 
outside the marked region) has not been affected at all by 
the filtering procedure. If one wishes to filter further the 
ground roll noise one may subtract successively higher 
modes. For example, in Fig. \7jp we show the filtered im- 
age after we also subtract the second eigenimage. One 
sees that there is some minor improvement, but remov- 
ing additional modes is not recommended for it starts to 
degrade relevant signal as well. 
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mode number 

FIG. 6: The relative energy of the first 10 KL modes of the 
region shown in Fig. 




FIG. 7: a) The filtered seismic section after removing the 
first eigenimage of the select region shown in Fig. |Sp,. In b) 
we show the result after removing the first two eigenimages. 



V. CONCLUSIONS 

An optimized filter based on the Karhunen-Loeve 
transform has been constructed for processing seismic 
data contaminated with coherent noise (ground roll). A 
great advantage of the KL filter lies in its local nature, 
meaning that only the contaminated region of the seismic 
record is processed by the filter, which allows the ground 
roll to be removed without distorting most of the reflec- 
tion signal. Another advantage is that it is an adaptative 
method in the sense the the signal is decomposed in an 
empirical basis obtained from the data itself. We have 
improved considerably the KL filter by introducing an 



optimization procedure whereby the ground roll region 
is selected so as to maximize an appropriately defined 
coherence index CI. We emphasize that our method, re- 
quire as input, only the generic alignment functions to be 
used in the optimization procedure as well as the number 
of eigenimages to be removed from the selected region. 
These may vary depending on the specific application at 
hand. However, once these choices are made, the filtering 
task can proceed in the computer in an automated way. 

Although our main motivation here has been suppress- 
ing coherent noise from seismic data, our method is by no 
means restricted to geophysical applications. In fact, we 
believe that the method may prove useful in other prob- 
lems in physics that require localizing coherent structures 
in an automated and more refined way. We are currently 
exploring further such possibilities. 
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APPENDIX A: RELATION BETWEEN THE KL 
TRANSFORM AND PROPER ORTHOGONAL 
DECOMPOSITION 

In dynamical systems the mathematical procedure 
akin to the KL transform is called the proper orthogonal 
decomposition (POD). In this context, one may view each 
column vector y*j of the data matrix A as a set of m mea- 
surements (real or numerical) of a given physical variable 
f(x,t) performed simultaneously at m space locations 
and at a certain time tj, that is, yjk = f(x = Xk, t = tj), 
k = 1, ...,m. For example, in turbulent flows the vectors 
iji often represent measurements of the fluid velocity at m 
points in space at a given time i. The data matrix A thus 
corresponds to an ensemble {j/j} of n such vectors, rep- 
resenting a sequence of m measurements over n instants 
of time. In POD one is usually concerned with finding 
a low-dimensional approximate description of the high- 
dimensional dynamical process at hand. This is done by 
finding an 'optimal' basis in which to expand (and then 
truncate) a typical vector y of the data ensemble. Such 
a basis is given by the eigenvectors of the time-averaged 
autocorrelation matrix R, which is proportional to the 
matrix T define above: 

1 n 1 

R^(yy t ) = -J2mv! = -r. (Al) 

»=l 

Hence the eigenvectors {Hi} of T are also eigenvectors 
of R. In POD parlance the eigenvectors {ui} are called 
empirical eigenvectors or proper orthogonal modes. In the 
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continuous case, the corresponding eigenfunctions Ui(x) 
of the autocorrelation operator are known as empirical 
orthogonal functions (EOF). 

From QJ, (JIJ and J5J), one can easily verify that 



Vr 



(A2) 



We thus see that the columns of the KL transform W 
correspond to the coordinates of the vectors y in the em- 
pirical basis: 



ki u k- 



(A3) 



k=l 



It is this expansion of any member of the ensemble in the 
empirical basis that is called the proper orthogonal de- 
composition or empirical orthogonal function expansion. 
It now follows from (IA3I) that 



(A4) 



x ' n ^ — ' n ^ — ' n ^ — ' 



where in the last equality we used the fact that 

= U*TU = A, (A5) 

where A is the diagonal matrix A = diag(Ai, A m ). 
Equation l|A4|l thus suggests that we can interpret the 
eigenvalue Xi as a measure of the energy in the ith em- 
pirical orthogonal mode. For example, in the case of tur- 
bulent flows where the vector yi contains velocity mea- 
surements at time i, the left hand of (|A4() yields twice 
the average kinetic energy per unit mass, so that |Aj 
gives the kinetic energy in the ith empirical orthogonal 
mode 0| . Similarly, in the case of seismic data the vec- 
tors yi represent amplitudes of the reflected waves, and 
hence the quantity J2i=i V? — J2i=i ma y be viewed as 
a measure of the total energy of the data, thus justifying 
the definition given in 

The optimality of the KL expansion also has a nice 
physical and geometrical interpretation, as follows. Sup- 
pose we write a vector y in an arbitrary orthonormal basis 



APPENDIX B: RELATION BETWEEN THE KL 
TRANSFORM AND PRINCIPAL COMPONENT 
ANALYSIS 

In statistical analysis of multivariate data, the KL 
transform is known as principal component analysis 
(PC A). In this case, one views the elements of a row 
vector Xi — (xn, Xi n ) of the data matrix A as being n 
realizations of a random variable Xi, so that the matrix 
A itself corresponds to n samples of a random vector X 
with to components: X — (Xi, X m ) , In other words, 
the column vectors ijj correspond to the samples of A. If 
the rows of A are centered, i.e., the variables Xi have zero 
mean, then the matrix T is proportional to the covariance 
matrix Sx of X [2(j : 

{Sx)ij = {XiXj) = ■ Xj = —Tij, (Bl) 



or alternatively in matrix notation 



Sx = (X X 



1. 



(B2) 



[Note that the matrices R and Sx defined respectively in 
(|A1|) and (|B2|) are essentially the same but have differ- 
ent interpretations.] In the PCA context, the diagonal 
elements Tu of the matrix T are thus proportional to the 
variance of the variables Xi, whereas the off-diagonal el- 
ements T^, i ^ j, are proportional to the covariance be- 
tween the variables Xi and Xj. Furthermore, the eigen- 
vectors Ui of r correspond to the principal axis of the 
covariance matrix Sx- The idea behind PCA is to intro- 
duce a new set of m variables Pi , each of which being a 
linear combination of the original variables Xi , such that 
these new variables are mutually uncorrelated. This is 
accomplished by projecting the vector X onto the prin- 
cipal directions of the covariance matrix. More precisely, 
we define the principal components Pi, i = l,...,m, by 
the following relation 



Pi = X ■ u t 



i=i 



HjXj. 



(B3) 



£« 



(A6) 



where Oj — e*i • y. If we now wish to approximate y by 
only its first k < to components, 



y k 



(A7) 



then the optimality of the KL expansion implies that 
the first k proper orthogonal modes capture more energy 
(on average) that the first k modes of any other basis. 
More precisely, the mean square distance (\y — y k \ 2 ) is 
minimum if we use the empirical basis. 



In other words, the vector of principal components P = 
(Pi, P m Y is obtained from a rotation of the original 
vector A: 



P = U l X. 



(B4) 



The covariance matrix Sp of the principal components is 
then given by 

S P = (PP^ = ^'il'f/j = -U l YU = -A, (B5) 
thus showing that 

{PiPj)=0, for (B6) 
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as desired. The first principal component P\ then rep- 
resents the particular linear combination of the origi- 
nal variables (among all possible such combinations 
that yield mutually uncorrelated variables) that has the 
largest variance, with the second principal component 
possessing the second largest variance, and so on. 

From J3J and l)B4j) one sees that the elements of the «th 



row of the KL transform W correspond to the n samples 
or scores of the ith principal component. That is, if we 
denote the sample vector of the ith principal component 
by Pi = (pa, —,Pin), then pij = • yj = For this 
reason in the PCA context the KL transform "3/ is called 
the matrix of scores. 
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